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A model system for classical fluids out of equilibrium, referred to as DPD solid (Dissipative 
fj , Particles Dynamics), is studied by analytical and simulation methods. The time evolution of a 

q_J ■ DPD particle is described by a fluctuating heat equation. This DPD solid with transport based 

on collisional transfer (high density mechanism) is complementary to the Lorentz gas with only 
kinetic transport (low density mechanism). Combination of both models covers the qualitative 
behavior of transport properties of classical fluids over the full density range. The heat diffusivity 
is calculated using a mean field theory, leading to a linear density dependence of this transport 
C/3 ' coefficient, which is exact at high densities. Subleading density corrections are obtained as well. At 

lower densities the model has a conductivity threshold below which heat conduction is absent. The 

e^ ' observed threshold is explained in terms of percolation diffusion on a random proximity network. 

The geometrical structure of this network is the same as in continuum percolation of completely 
I , overlapping spheres, but the dynamics on this network differs from continuum percolation diffusion. 

' O ' Furthermore, the kinetic theory for DPD is extended to the generalized hydrodynamic regime, 

C ' where the wave number dependent decay rates of the Fourier modes of the energy and temperature 

O fields are calculated. 
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05. 40. -a Fluctuation phenomena, random processes, noise, and Brownian motion 
64.60.Ak renormaUzation-group, fractal and percolation studies of phase transitions 
05. 10. -a Computational methods in statistical physics and non-linear dynamics 



I. INTRODUCTION 



(N 
> 

(N 

O 

\o 
o 
^. 

^p ■ The Lorentz gas, describing the self diffusion of a moving particle in a random array of scatterers, 

has played an important role in understanding the transport properties of classical fluids, and in 
developing and quantifying the role of correlated sequences of binary colhsions (ring collisions), as 
well as in extending the kinetic theory to moderately dense fluids (see [H H, |M H, Q and references 
™|j . there in). 

r^ ' The Lorentz gas contains only the mechanism of kinetic transport, which is the most important 

Q , transport mechanism at low densities. In dense fluids there is a second mechanism for transport, 

O ' called collisional transfer ^ , through which energy and momentum is instantaneously transferred 

through the interactions between particles within each others force range. This is the dominant 
mechanism at high densities. 

In this paper we discuss a model, complementary to the Lorentz gas, which contains only the 
mechanism of collisional transfer, and for which transport properties can be evaluated exactly in 
the high density limit. It is also important to develop systematic theories for subleading large 
density corrections. The combination of both complementary models, the Lorentz gas and the 
DPD solid, might be able to describe the qualitative density dependence of transport properties of 
classical fluids over the full range from low to high densities. 

Before introducing the random DPD solid we briefly discuss the relevance of DPD models for 
the study of classical fluids. This stems from the great interest in computer simulations of complex 
fluids and quenched random media, which are challenging problems as several different space and 
time scales may be involved. Fully atomistic simulations of such systems fail in reaching the 
larger scales, and different mesoscopic models/techniques, such as Dissipative Particle Dynamics 
(DPD), Smooth Particle Hydrodynamics, cellular automata and lattice gases, lattice Boltzmann 
methods, etc, offer possibilities to explore these larger scales. For that reason the DPD technique 
was originally introduced Q as a mesoscopic particle method for simulating complex fluids and 
colloidal suspensions. It is therefore important to provide a theoretical analysis of such systems, 
as will be done in this paper. 

The idea of the method is that each DPD (ooint) oarticle reoresents a mesoscooic Dortion of 



fluid. The interactions among these point particles have no hard core, and are softly repulsive. The 
lack of hard core interactions allows time-step driven algorithms 8] . In the original formulation 
the DPD particle is described in terms of its position and velocity with three different types of 
interactions: conservative, dissipative and stochastic. The forces between particles are pair- wise, 
such that mass and momentum are conserved, but energy is not. This formulation is restricted to 
isothermal problems, but describes a proper hydrodynamic behavior for viscous fluids 8, ft llfllllj l 
in a large number of pro blems. The method has also been successful in describing properties of 
colloidal suspensions [Ij, U^ , polymer solutions [ij, ^^ , phase separation [l^ [13 or membranes 
[l8lll9l |. In standard DPD the forces are Gaussian white noise, which have been recently extended 
to colored noise 20]. 

A generalization of the model to include energy conservation has been proposed as well [22, l2a | 
in order to describe heat flows and other thermal effects in fluids out of equilibrium. In the picture 
where DPD particles are understood as droplets or mesoscopic clusters of microscopic particles, 
one can consider the kinetic energy lost in dissipative interactions as being transformed into energy 
of internal degrees of freedom of a particle. The number of internal states of a DPD particle with 
energy e, exp[s(e)/fcs], is modelled by an entropy function s(e), which implies a temperature T{e) 
defined through ds{e)/de — 1/T{e). The evolution of the internal energy has two contributions. 
The first one describes how the friction forces contribute to the change of kinetic energy; this 
is viscous heating. In addition, the phenomenon of heat conduction has to be considered, where 
energy or temperature differences between particles (subsystems) produce a flux of internal energy. 
In Ref. 23] Bonet and Mackie have estimated the transport coefficients in the limit of high friction 
for this thermal DPD fluid, using a method somewhat similar to that used in Ref. 'flj- In R-^f. 
|24| we have investigated the DPD fluid obeying the conservation laws of mass, momentum and 
energy, and we have calculated the full set of transport coefficients in the Navier-Stokes and energy 
balance equation, including kinetic and coUisional transfer contributions, as well as their wave 
length dependent generalizations. However, as pointed out in Ref.[3 for standard DPD fluids, the 
theoretical values for the transport coefficients appear to agree only well with computer simulations 
at higher densities. 

To analyze the difficulties and to develop a more accurate description it is of interest to consider 
a simplified model with heat conduction, the random DPD solid, which still has the basic features 
of DPD. Such a model can be obtain by considering the energy conserving DPD models |2lll2^. 
where the dynamical degrees of freedom of a particle, Xi{t) — {ri{t), Vi(t), ei{t)} (i = 1, 2, • • • , A'') 
are position, velocity and internal energy. By freezing/quenching the velocities, the particles can 
be characterized by static (random) positions r.^, and by dynamic energy variables ei(i), with their 
total energy, E — ^- ei{t), conserved. 

Consequently, the macroscopic evolution equation for this system is Fourier's law of heat diffu- 
sion, and the system is able to carry a macroscopic heat flux, provided a macroscopic fraction of 
the particles is inside each other's interaction ranges. For simplicity we set the conservative forces 
equal to zero (point particles), and take the interaction range of the dissipative and stochastic 
forces equal to re- This model has been proposed in J25j ]. basically as a discrete fluctuating heat 
equation. 

Transport of energy in fluids consists in general of kinetic transport, carried by moving particles, 
and instantaneous transport through the interactions, the so-called coUisional transfer. In the 
present model there is no kinetic transport, and the only type of transport is through coUisional 
transfer. In dense fluids coUisional transfer is also the dominant mechanism of transport. 

The basic observation is that the coUisional transfer mechanism can be viewed as hopping of 
energy across existing bonds, as illustrated directly by the dissipative part of the A^— particle 
Langevin equations, i.e., 

dei/dt = \o2_^w{rij){ei- ej). (I.l) 

3 

It is essentially a discrete diffusion equation on a random network with bonds, to be deflned below. 
In the previous equation Aq is a relaxation coefficient, and w{rij) is a positive interaction function, 
say, equal to 1 for rij < re and elsewhere, where re is the interaction range. Here we are dealing 
with a dynamic diffusion problem on an underlying static percolating structure. In order to have 
any transport of energy the density should be sufficiently large, such that there exists a percolating 
path of connected particles, spanning two opposite boundaries of the system. Groups of isolated 
or connected oarticles. which are not oart of the oercolatinff oath, form non-conductine islands. 



So the underlying static structure is a bond-percolation cluster, whose density should not be below 
some threshold density |26l |. 

The question is then, what is a bond in this quenched random solid of point particles occupy- 
ing random positions {vi\i = 1,2,- • • , iV} and having an interaction range re- To visualize the 
connected network, we connect every pair {i,j} of fixed point particles with a bond if \rij\ < re- 
Energy can only hop between connected particles. With this definition of a bond, we have a well 
defined percolation diffusion model on a random proximity network with a constant hopping rate 
per bond. 

The geometrical structure of this connected network is the same as in continuum percolation 
of completely overlapping spheres ( see Ref. [271 l2q and references therein). Suppose we put 
black circles of radius R = \rc on every point particle. Then two particles {i,j} are considered to 
be connected or "overlapping" if rij < 2R = re, and we obtain the above continuum percolating 
structure. 

However the dynamics of the present diffusion problem is quite different from continuum per- 
colation diffusion models, such as the overlapping Lorentz gas Q and the Swiss Cheese model, 
where diffusion occurs in the void spaces outside the overlapping spheres [23, [SQ] , or the Inverted 
Swiss Cheese model and others 29] , where diffusion occurs in the complementary space, i.e. inside 
the overlapping spheres. Such models can be mapped either on the discrete random proximity 
network (Inverted Swiss Cheese model), or on its dual network (overlapping Lorentz gas, Swiss 
Cheese model). The big difference is that the maps of the continuum diffusion models have a 
wide distribution of hopping rates, usually singular at small rates, because of the appearance of 
bottleneck passages [Slj , whereas the random DPD solid has constant or nearly constant hopping 
rates (depending on the shape of the range function w{r)). 

We also note that the overlapping Lorentz gas and the Swiss Cheese model are percolating below 
a threshold density, whereas the dual models: the Inverted Swiss Cheese model and the random 
DPD solid are percolating above a threshold density, illustrating the relevance of these models for 
low or high density fluids. 

In the overlapping Lorentz gas, say in two-dimensions, the disordered network is formed by the 
vertices and edges ("bonds") of the polygons that partition space into a Voronoi tessellation [5ll29l|. 
Here the blocked bonds are the edges that are perpendicular to and bisect the lines of centers with a 
length rij < re- Hence, the blocked bonds in the Voronoi tessellation are the duals of the bonds in 
the random proximity network. Dynamical properties (exponents, amplitudes) near the threshold 
density will in general be different on discrete disordered networks with constant hopping rates, 
such as the random DPD solid, and on continuum percolation models, corresponding to networks 
with a wide (singular) distribution of hopping rates. 

From the point of view of dense fluids the percolation phenomena are in a way just an interesting 
pathology of the random DPD solid, caused by freezing out the translational degrees of freedom 
of the corresponding DPD fluid. So the main interest of this paper is focused on a quantitative 
description of transport coefficients away from the percolation density. This kinetic theory problem 
has not received much attention in the literature of the last three decades, which has been focusing 
on the behavior near the threshold, and not on the density dependence away from the threshold 
density. 

So far we have described the dissipative part of the fluctuating heat equation. An A^— particle 
state described by the dissipative equation (|I.1|I would decay from an arbitrary initial state to a state 
of zero temperature with all energies e^ — E/N. To make the DPD solid reach thermodynamic 
equilibrium one adds a fluctuating force to the evolution equation, satisfying the fluctuation- 
dissipation theorem. How this is done is explained in Section 2 and Appendix A. 

The plan of the paper is as follows. It starts with a more detailed presentation of the heat 
conducting random DPD solid in Section II, while some more general properties, such as the Ti. 
theorem and the equilibrium properties are discussed in Appendix A. In Section III A the heat 
diffusivity is calculated using a mean field theory which is expected to be exact at large densities, 
where local density fluctuations can be neglected [ifl HjI ■ In Section III B the large deviations 
at low densities between the results of mean field theory and simulation are explained in terms 
of bond percolation on a random proximity network. In Section IV we derive the wave number 
dependent decay rates of the spatial Fourier modes of the fluctuations in energy and temperature 
fields. The last section is devoted to some conclusions. 



II. DPD - HEAT CONDUCTION MODEL 

The heat conduction in dissipative particle dynamics is modelled as a thermally conducting 
quenched random solid. The system is described by TV = nV point particles at quenched random 
positions r^, contained in a volume V — L"^. Each DPD particle is a mesoscopic subsystem, that 
interacts with all particles that are inside its interaction sphere of radius re- The only dynamical 
variable is the internal energy of the particle, e^, which captures the internal deg rees of freedom of 
the mesoscopic particle. Its evolution equation is the Langevin equation |25ll3^ . 



de, - ^'Afo') (Tj - T,) dt + J2' '^iiJ)dW,j (t), (II.l) 



where the prime indicates the constraint j ^ i. The first term on the r.h.s. is dissipative and 
specifies that a temperature difference causes flow of energy. The second term represents the 
Langevin force, described as Gaussian white noise, 

Fij {t)dt = a{ij)dWi.i (*)• (II-2) 

It takes into account thermally induced fluctuations in each particle causing random exchange of 
energy between particles. Conservative forces are absent in this model. The relaxation coefficient 
X(ij) models thermal conduction and a(ij) is the amplitude of the noise. These model parameters 
are assumed to be symmetric under particle interchange. In principle X{ij) and a{ij) depend on 
the relative distance r.^ , and on the internal energy of the particles i and j. If a{ij) depends on e^ 
and Cj, then HII.2() represents multiplicative noise, because the internal energies themselves depend 
also on the noise. The model parameters are of the general form, 

X{ij) = Xtjw{rij); a{ij) = a^jVUsirij). (II. 3) 

The range or interaction functions w(r) are Ws{r) are positive, they vanish for r > re and have a 
finite non- vanishing value at the origin. Moreover, we choose the following normalization for w(r). 



w] = / drw{r) = rf (II.4) 



and we will see below that the relation Ws(r) — y'wlr) is imposed by the detailed balance condi- 
tions. 

To define the temperature Ti of a mesoscopic particle with energy e^ the equation of state has 
to be specified, for which we use the entropy, or equivalently, the density of internal states. The 
simplest choice is the entropy for an ideal solid, 

s(e) = ain(e/e„), (II.5) 

where Cy = aks is the heat capacity of a DPD particle, which is assumed to be a constant, 
independent of e. The parameter e^ is a constant reference energy, and represents an additive 
constant to the entropy. Consequently, the temperature of a DPD particle follows from ds/de — 
l/r(e), and is given by, 

ei = CyT{ei) ^ aksTi. (II.6) 

The dimensionless number a = Cy/ks is a measure of the size of the DPD particle because it 
scales like the number of internal degrees of freedom of the particle. Hence, a is a large number, 
and subleading corrections of relative order 1/a will be consistently neglected in this paper. 

The total energy, E = X^i^*' ^^^ only internal energy contributions, and is exactly conserved 
by construction, because the r.h.s. of ljll.l|) . when summed over «, is chosen to be anti-symmetric 
under particle interchange, and therefore vanishes. Consequently, the increments of the Wiener 
process associated with the heat conduction have to be antisymmetric dWij = —dWji. The noise 
term represents Gaussian white noise with a mean dW = and with a noise strength. 



dW,At)dW.u-(s) = (6,,>S,.- - S,.'S,,>) mi-aidt.ds\. (11.7) 



where dWdW represents an average over the random noise, and minja, 6} denotes the minimum 
of a and b. The corresponding Fokker Planck equation for the time evolution of the N particle 
distribution, p(X, t) in the phase space given by X = {x, = (r^, ei)\i = 1,2 ■■ ■ N}, reads 



dtp = Lp. 



(11. 



If the stochastic differential is interpreted according to Ito |34L ISSL ISo. ISTj the Fokker Planck or 
Liouville operator L and its adjoint U have to be written as, 



L = J29k 
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\izj){T,-T,) + -d,,a^zj) 



\{i]){T,-T{) + -a\ij)d,, 



^^i, 



(11.9) 



where dij = d/dci — d/dej. Consequently the Fokker Planck operator has the standard form, 
appropriate for a multi-particle Fokker Planck equation. In Appendix A we construct an Ti. function 
and prove an Ti. theorem, i.e. dtTi. < 0, which holds provided the dissipation coefficient X(ij) and 
the noise strength a(ij) satisfy the so-called detailed balance conditions, 

a'{i])^2kBm)T,T, and d,ja^{ij)^Q (V{i,j}), (11.10) 

which implies that w1{r) — w{r), as derived in Appendix A. Inserting this relation into 1)11.1(1 puts 
the Langevin equation into the form. 



dT,,^ 
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■dci = 



-^y' m){Tj - T,)dt + J2kBX{ij)T,T, dW,, 

OiKB "^"^ ' 



(11.11) 



where the fluctuating term represents multiplicative noise. To further specify \ij we consider 
a subsystem 1 (here a mesoscopic particle) with energy ei = CyTi, in contact with a second 
subsystem of temperature T2. Then irreversible thermodynamics gives for the energy relaxation 
in subsystem 1, 



dt [cj dt akB^ ^ 



Ti 



(11.12) 



where the specific heat of the system is C^, — akg, and the relaxation coefficient A12 — akgXQ 
is a material constant, independent of the energies {ei,ej} of the interacting subsystems and 
proportional to the size a of subsystem 1. So, we choose 



A,: 



aksXa, 



(11.13) 



which defines the relaxation parameter in the Langevin equation for our heat conduction model, 
which reads finally. 



dTi = ^ Xow{rij){Tj - Ti)dt + y 2a ^Xaw{rij)TiTj dWi^ 



(11.14) 



where Ti, Tj may be replaced by e.;, ej. When performing simulations it is convenient to make the 
variables in the Langevin equation dimensionless, i.e. we express distances as r = r/vc, the time 
as i = XqI, the increment of the Wiener process as dWij = \/XodWij, and the temperature as 
Ti — Ti/Tu, where T„ is an arbitrary reference temperature. It follows by setting Aq = 1 in the 
equation above. Here dWij satisfies the relation (|II.7p with dt and ds replaced by dt and ds 
The corresponding Fokker Planck operators take the form. 



L = afcsAo ^ w{rij)di.j [Ti - Tj + ksd^jT^Tj 



t<] 



L^ = aks Xoy2 w{rij)[Tj -Ti + kBTiTjdijjdi;^ 



(IL15) 



From the discussion in Appendix A leading to (|A.10|I it follows that a?{ij) and dij in ljll.91) . 
and TiTj and dij in (jII.15|l are commuting to leading order in 1/a. This implies that the Ito and 
Stratonovich interpretations are coinciding, and that the stochastic differential dWij can be treated 
as a differential in ordinary differential calculus. 

In the present context the Fokker-Planck equation is frequently referred to as Liouville equation. 
In the same spirit the corresponding mesoscopic Langevin equations are referred to as microscopic 
equations. To complete the microscopic description we derive the local conservation law for the 
microscopic energy density. This supplies us with a microscopic expression for the energy flux. It 
will be used in the simulations to measure the macroscopic heat current and the heat conductivity. 

We first introduce the static (quenched) particle density and the dynamic energy density, given 
respectively by, 

n(r)=^,5(r-r,); e(r) = ^e,J(r ~ r,). (11.16) 

i i 

Following standard arguments we derive an expression for the local m,icroscopic energy flux g(r) 
as well as for the total flux Q — J drq{v) . A hat on a symbol denotes a mesoscopic quantity. 
The equation of motion for the expectation value (e(r))t — e(r, t) becomes, 

dte{r,t) = J dXe{r\X)dtp = {L^e{v))t = -V • (q(r)),. (11.17) 

The last equality in ljII.17|) shows the local energy conservation law. We have further used the 
relation, 

^^e(r) = Y.^<, w{n,)X,,{T, ~ T{) [5{v - r,) - 5{v - r,)] 

^ -V • E.<, w{n,)v,,K,{T, ~ Ti)5{v - r,) ^ -V • q(r) , (11.18) 



and note that the terms in the Fokker Planck operator (|II.15|) . coming fronr the noise, do not 
contribute to the macroscopic flux of energy. The last line has been obtained by expanding 5{y~Tj) 
in powers of r^ , i.e. 5{y — r^ + r^) ~ 5{y — r^) + Ty • V5(r — r^) + O(V^). The total microscopic 

heat flux, Q = / rfrq(r) = Qjj + Q^, consists of a dissipative (D) and a fluctuating (R) part |3q . 
With the help of the relations Ay = afc^Ao and e^ — aksTi the dissipative part is given by, 

Qr> = Ao ^ w{rij)rij{ej - e,,), (11.19) 

where terms of order 1/a have been neglected. We note that the current Qd does not contain 
kinetic contributions, but is a sum of pair contributions involving the dissipative interactions. This 
is the mechanism of coUisional transfer, representing instantaneous transfer of energy to particle 
i from all particles j in the interaction sphere defined by Vij < Vc- The current Qu should 
be compared with the contributions to the microscopic stress tensor involving the conservative 
interaction potentials in Hamiltonian fluids. Furthermore, the expression for Qd also illustrates 
that the total macroscopic heat current is determined by the energy difference, ej — e^, i.e. by the 
"temperature gradient" between r^ and Tj. The fluctuating part Q^, = — X]i<7 «(*i)i"jj-f'(u) with 

Fij in Eq. lITOl and (fiR)t = 0. 

The macroscopic energy flux q obeys the standard linear constitutive law of irreversible thermo- 
dynamics, 

q=(q)t = -AVT, (11.20) 

where A is the coefficient of the heat conductivity. Combination of (|II.17ll and ljll.20|) with the 
relation 6e = nCyST = naksST, yields the equation for heat diffusion, 

dtT = — V^T EE DtV^T, (11.21) 

naks 

where Dt is the heat diffusivity. 

For later reference we mention that the heat conductivity can also be calculated and/or simulated 



akBn{r)T{r,t), satisfies a local conservation equation, the heat conductivity can be expressed in 
the Green-Kubo formula, 

where the average (• • • )o is taken over the apropriate equilibrium ensemble. On the basis of the 
analogy between the Liouville and the Fokker Planck operators in IjII.lSp . with Tj replaced by 
ei/aks, any of the standard derivations of these formulas |33,ll3| carries over directly to our DPD 
solid. We further note that the microscopic energy current Q(t) does not contain any "subtracted 
part" because this model does not have any conserved quantity with a vector character, such as 
the total momentum. 

For the case of general DPD fluids the Green-Kubo formula for the viscosity has been derived 
in Rcf. [JJ. Another representation of the transport coefficients, equivalent to the Green-Kubo 
formulas, is given by the so-called Helfand formulas |4^i4^. It reads for the present case, 

with M given by, 



One easily verifies that the microscopic heat current, Q^ — L^M, in ljll.22(l can be obtained from M. 
The Helfand formulas are generalizations of Einstein's formula for the coefficient of self diffusion, 
and are presumably more convenient in numerical simulations than the Green-Kubo formulas. 

A further consequence of the Ti. theorem, discussed in Appendix A, is the existence of a unique 
equilibrium state, the Gibbs' state. Its explicit form has also been determined in Appendix A. In 
the main text of this paper we only need the single particle equilibrium distribution function for 
the DPD sohd, as derived in HA.15|) . i.e. 

^o(^) = r( In (^^)" exp[-/3e], (11.25) 

T[a + 1) 



where a is a measure for the number of internal degrees of freedom. In Ref. [23, I24L I2!t| simulations 
of the equilibrium distributions in the conduction model have been performed, and good agreement 
between the simulation results and the analytical expressions has been obtained. For instance, the 
simulated energy fluctuations agree very well with the theoretical prediction {{5tiY) = (a-|-l)//3^ ~ 
a/ 0^ within error bars smaller the 0.7% |2^. 



III. HEAT CONDUCTIVITY 
A. Mean field theory 

The DPD model for heat conduction is expected to produce a macroscopic behavior described 
by a macroscopic heat equation. Our aim is to prove this assertion and to relate the effective 
thermal diffusivity appearing in the macroscopic heat equation to the model parameter Aq and the 
range function w{r). To do so we will use a mean field approximation. We start with an a priori 
estimate of the transport coefficient. 

In a naive kinetic picture of the relevant transport mechanism, used in Ref. [23 , amounts of heat 
or energy hop on a random lattice with an average lattice distance Is = n~^''^, and a hopping 
frequency wq- This picture, based on kinetic transport of energy, leads to a heat diffusivity Dq = 
wo's J where wq is the decay rate of an energy or temperature fluctuation. As the decay rate loq (x n 
(see below), this would lead to an a priori estimate Dq cx n^~'^/-^, which does not hold for the 
collisional transfer mechanism. 

As the DPD particles are quenched, there is no kinetic transport, but only collisional transfer of 
energy, i.e. instantaneous transfer of energy through particle interactions. It takes place only over 



of 13 = ujqt^, where Tc is the range of the interaction function w{R), and wq is a typical frequency. 
For large densities this frequency can be estimated from the first term on the r.h.s. of ljII.14p as, 



Wo 



AoX! (^(^y)) - ^on[w] = pAo, (III.l) 



where [w] is defined in (|II.4p and (• • ■ ) denotes an average over all quenched particles. Here the 
reduced density p = nrf is on the order of the mean number of interacting neighbors with r^ < rc^ 
surrounding the i— th particle. The freedom to shift constant factors in (|II.3|I from Ag to w{r) is 
fixed by the normalization \w\ — rf in ljll.4|) . The a priori estimate for the diffusivity at large 
densities is then, 

D-^y^^pxy^. (III.2) 



The estimate (jIII.2|l will be confirmed by detailed kinetic theory calculations. 

To calculate the heat flux we express the average of (jII.19|l in terms of the two-particle distri- 
bution function, /(^^(ri, ei, r2, e2, i), yielding 

q - ^Ao /" deidez / dRu;(i?)R(e2 - ei)/*'Hr, ei, r - R, €3, t). (III.3) 

The basic ansatz in this mean field theory is that the fluid rapidly relaxes to a local equilibrium de- 
scribed by the local fields 6(r, t), being the temperature T(r, t) = l/fcB/3(r, t) and local (quenched) 
density n(r). This happens on a time scale I/loq where ujq ~ Aop(r) is the local relaxation rate of 
the temperature fluctuations and p(r) is on the order of the number of particles inside the sphere 
centered at r. Temperature gradients, which are smoothly varying in space, only build up on 
a hydrodynamic time scale, as described by Fourier's heat law. Because conservative forces are 
absent in our model, the local equilibrium pair distribution function is exactly equal to the corre- 
sponding pair function of an ideal gas, i.e. it is simply a product of single particle local equilibrium 
distribution functions /o, i.e. 

/(2)(r, e, r', e\ t) = fo{e\b{r, t))/o(e'|5(r', t)), (IIL4) 



which depends on 6(r,t) and b{r',t). Their explicit form follows from ljll.25|l as, 

/o(e|6(r)) = n{r)M^\f3{r)) = n(r)/?(r) [/3(r)e]" e-^W7r(a + 1), (III.5) 

where ipo only depends on the temperature at the position r of the particle. 



To derive an expression for the heat conductivity in ljll.20(l the equations ljIII.3|l - (|III.5(l need 
to be expanded in gradients, i.e. /3(r — R) = /? — R • V/9 + ©(V^), and similarly for n(r — R). 
Because the integrand in (|III.3p is anti- symmetric in both R and (e2 — ei), the only non- vanishing 
contribution to q from /^^^ in (|111.4() must be linear in R and 62. The result is, 

/o(e2|/3(r - R)) ^ -n(R • VP)M^2)-^ In ^0(62) ^ ne2(R ■ V/3)Vo(e2). (III.6) 

Symmetrizing the last expression, £2 — * ^(^2 ~ ei), and inserting this in ljIII.4p and (jIII.Sp yields 
after some algebra for the heat current, 

q = -iAofcsn^ / dRu>(i?)RR • VT /3'^{{e2 - ei)2)o = -AVT, (III.7) 

yielding a mean field prediction for the heat conductivity, 

A(p) = ^XokBn^[w]{R^)^ = Xooip), (III.8) 



which is quadratic in the density. We refer to this mean field value IjIII.Sp as the saturation or high 
density value. This value is expected to be exact when the difference between the actual number of 
particles inside an interaction sphere and its mean value can be neglected. Here ((£2— ei)^) ~ 2q;//3^ 
denotes an average over the canonical ensemble (|A.11|I . and is given by (|A.16I) . Moreover the 
following relation has been used. 



I ^VtR P^i/./RA — _A „ rlTiJi^^n(m = -/i ^U„l/f?2\ 



riTT Q\ 



where the last equahty defines the normahzed second moment {R'^).^. The mean field value for the 
heat diffusivity, defined in (jll.21|l . is obtained similarly, 

Doo = ^^^{R')n., (III.10) 

where uq = Ao«.[w] = Xqp is the typical decay rate of an energy fluctuation. 

In the literature on DPD simulations different choices of the range function w{r) have been 
considered, for instance, 

(Heaviside) 

u,{r) = { ^n' " ~y "^''^ ~ ''' „ (Triangle) (jjj_^^) 

; - r) (Lucy) 

Here 0{x) is the Heaviside or unit step function, and the constants Ad, Bd, Cd are normalized such 
that [w] ~ rf in d— dimensions. Calculation of the moments {R^)w yields then, 

2(d+2) ^o?^c (Heaviside) 

^- = 3^ - { 2(d+t)(rf+3) ^org (Triangle) (HI.12) 
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Before concluding this subsection we make a number of comments. If we would have chosen 
a normalization of w(r), different from llll.4|) . say w'{r) = Cw{r), while keeping the unit of time 
fixed, the diffusivity would change to D' = CD. We recall that the heat current in (jIII.311 . which 
is based on the mechanism of collisional transfer, establishes itself instantaneously, i.e. a high 
frequency limit occurring on the fast time scale l/^Oj because pair interactions transfer energy 
instantaneously between r^ and r^, whenever r^ < re- Such heat currents only occur in dense 
systems, and they are non-vanishing in a state of local equilibrium. On the other hand kinetic 
transport of energy, carried by moving particles, establishes itself only on the slower hydrodynamic 
time scale. 

We add a technical remark. In general the heat flux in (|III.3|) would also pick up contributions 
from the additional gradient terms in the single particle distribution functions in ljIII.4|) . 

fix^) = foix^)[l + iJ(e,|/?(r,))V/3(r,) + •••]. (IIL13) 



with i = 1,2, that would have to be added to lllll.4|l . However, such terms contribute to the heat 

flux q only terms of ©((VT)^) and O(V^T) because of the parity in R of the integrand in ljIII.3|) . 

Next we note that the results (|III.10|) and HIII.12() above are in qualitative agreement with the a 

priori estimate, D ~ wor^ in l|III.2() . The expressions for Dt in IjIILlOjl and (|III.12|I . with physical 

dimensions L'^ /t, has the typical structure of a diffusivity, i.e. a collision frequency multiplied 

with the square of the interaction range re, over which energy is transported by the mechanism of 

collisional transfer. 

Finally we note that the heat diffusivity and the typical frequency wo are proportional to the 

reduced density, p = nrf, whereas the heat conductivity Aoo in (|HL8|I is proportional to p^. In the 

DPD fluid with viscous dissipation, studied in Ref. [lOj, there are of course kinetic contributions to 

the viscosity as well, apart from the collisional transfer (ct) contribution to the kinematic viscosity, 

^ct ~ P, being a diffusivity, and to the shear viscosity, 77^4 ~ pvct ~ P^- Also in the Enskog theory 

for a dense fluid of hard spheres |g] similar collisional transfer contributions to the transport 

coefficients are present, sometimes referred to as the instantaneous transport coefficients, 7700, Aoo, 

2 



for the reasons explained above. The p — density dependence of the heat conductivity in ljIII.8|l is 
a direct consequence of the collisional transfer mechanism. For sufficiently high densities, the heat 
flux q is proportional to the local density of interacting pairs, ~ p^, and the heat conductivity 
is given by its saturation value Aoo(p) in (|HL8|I . This prediction holds when the typical density 
fluctuations Sp ^ ^/p can be neglected with respect to the mean value p. 

B. Simulations and conductivity threshold 



In our simulations we use the Langevin equation ljII.14|) to determine the dynamic properties of 
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with a cold and a hot waU at x = and x — L respectively. This is achieved by putting two extra 
layers of particles at each boundary, filled with particles at a constant temperature. The layers have 
a width Tc to ensure that any particle inside the system interacts on average with the same number 
of neighboring particles. Therefore, these extra layers act as thermal baths which are prepared at 
temperatures Tcom and Thot such that a gradient (Thot — Tcoid)/L is applied to the system. In the 
remaining directions we impose periodic boundary conditions. Initially the box is seeded with N 
mesoscopic (point) particles located at random, surrounded by overlapping interaction spheres of 
radius re- The initial temperature of the particles is (Thot — T'coid)/2. Another relevant quantity is 
the mean num,her of interacting neighbors^ ^{'"'c) — (yi) inside the sphere rij < re, where 



l^i = PiWd 



E'^(^ 



iy{rc) ~ nujdrf = prud (large p). (III. 14) 

Here p is the reduced density, and zud — 7r'*/^/r(l + d/2) is the volume of a d— dimensional unit 
sphere {d — 1, 2, • • • ). The fluctuating variables Vi = vjdPi are subject to large static fluctuations, 
especially at small densities, and are distributed according to a Poisson distribution. 

We perform a series of simulations by numerically integrating the equation (|II.14p with a given 
temperature gradient and compute the macroscopic heat flux in the steady state. The time re- 
quired to reach the steady state increases with decreasing density, and even diverges as the density 
approaches a threshold value, to be discussed below. The heat flux is calculated with the help of 
deflnition ljII.19|l . This expression involves a large number of pairs of particles which guarantees 
reasonably good statistics. If the density is not too close to the threshold, the simulations show a 
linear temperature profile between the two heat baths, as predicted by Fourier's law. The measured 
heat flux is found to be linear in the applied temperature gradient, from which the heat conductiv- 
ity can be extracted. Therefore, since this linear relation between Qx and V^jT has been conflrmed 
by the simulations, a single measurement of the macroscopic heat flux for a given gradient suffices 
to obtain the thermal diffusivity. Statistical errors can be estimated by making an average over 
several independent runs. In principle, alternative ways to measure the heat conductivity would 
be to simulate the Green-Kubo formula (|II.22p . or to set up a sinusoidal intitial temperature field 
Ti(0), and to measure Dt from the decay rate of the initial temperatures p4| . 

For our further discussions it is convenient to introduce, 

Rd{p) = Dt{p)/D^ = A(p)/Aoo(p). (III.15) 

We first consider the transport properties in the three-dimensional heat conduction model, and use 
the simulation results, obtained in Ref. [2J,|23, where the range function w{r) was chosen to be the 
Lucy function Hlll.ll() , and we compare the measurements with the analytic predictions H111.12() , as 
shown in Fig.l. The simulation results were performed with random spatial configurations. They 
approach the theoretical predictions at high densities. However, as the density decreases, the ratio 
i?3 rapidly decreases, almost by a factor 2 at the lowest densities simulated {p ~ 3.8). In these 
measurements averages over different values of the parameters N and a have been taken at fixed 
reduced density p. Both N and a should be sufficiently large to reduce the finite size effects and 
improve the statistics. 

A consistent explanation for these large deviations seems to be that the mean field theory does 
not take into account the local fluctuations in the actual number of particles, pi, in the interaction 
sphere around particle r^. These fluctuations are particularly large at low densities. To test 
this working hypothesis we place all N particles on a completely filled cubic lattice with lattice 
distance Ig ~ {V/NY^^ , thus suppressing all local density fluctuations. The resulting measurements 
are represented in Fig.l by (A)'s. Note that this suppression of density fiuctuations considerably 
extends the agreement between theory and simulations towards lower densities. The on-lattice 
simulations support our working hypothesis, and the observations are consistent with the good 
agreement at high densities, where fiuctuations in pi are small, but a theoretical explanation of the 
p-dependence of the heat conductivity of our original random solid is still lacking. The improved 
agreement between theory and simulations was here obtained by modifying the model. Further 
modifications of the random heat conducting solid model to suppress the local density fluctuations 
were introduced in Ref. |45j | , but do not increase our understanding of the density dependence of 
Dt{p) in DPD fluids and solids. 

To test these concepts the simulation results for the three-dimensional model in Fig.l are not 
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Figure 1: Simulated values of the thermal diffusivity Rsip) — Dt/ Dao versus p for the 3-D heat conduction 
model with A'' = 10^ DPD particles. Results from off-lattice simulations (•) and on-lattice simulations 
(A) are compared with the mean field value (IIII.121 (dashed line). At p = 3 and 27 the system sizes are 
respectively L/tc ~ [N/py''^ ~ 6.9 and 3.3, which is rather small, and strong finite size effects are to be 
expected. 
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Figure 2: 2-D simulations of the dimensionless heat diffusivity DT{p)/\or^ plotted versus p to test the 



mean field prediction Doo{p) in 1111.811 (dashed line), valid at high density. Labels (A) and (•) refer to the 
random solid respectively with and without Langevin force in systems with 10^ and 10^ DPD particles. H 
refers to the Heaviside weight function. 



(L/tc — 3.3 or 6.9). Furthermore the range function w(r) was taken to be the Lucy function 
in Hlll.ll() . which gives larger weights to shorter bonds. To optimize the similarity with the 
classic bond percolation problems, we give all bonds equal weights by taking w(r) as the Heaviside 
function. Moreover, to make the simulations less demanding, we have carried out simulations 
of our heat conduction model in two dimensions. Fig. 2 shows the good agreement between the 
two-dimensional simulations and mean field theory at high densities, as also observed in the three- 
dimensional simulations both of Fig.l, as well as in Fig.l of Ref. |4g. However, at lower densities 
Dt{p) decreases faster than linear, as shown in Fig. 3. 

To display the connection to percolation it is instructive to plot i?2(p) = DT{p)/Dao{p), as 
shown in Fig. 3. The plots strongly suggest the existence of a conductivity threshold pc > Pp, 
where pp = 1.43629 is the threshold value of two-dimensional continuum percolation |27l |. 

Fig. 3 shows that in our model the conductivity threshold 1.4 < pc < 1.5 which agrees quite 
well with the known percolation threshold. Nevertheless, determining the value of pc with higher 
precision becomes more delicate since the required times for equilibrating the system are diverging 
as p J, Pc- Simulations show that the approach to the expected linear temperature profiles at p > 4 
is relatively fast (relaxation time tg < 300). But these times increase for smaller densities, for 
p = 1.6 {to — 2.5 X 10^) and for p = 1.4 {to — 10^). The profiles at the times of measurement are 
shown in Fig. 4. Furthermore Figs. 3b and 3c show the very slow cross-over of the conductivity at 
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Figure 3: 2-D simulations of the heat diffusivity R2{p) plotted versus p. Figures a,b,c show respectively the 
behavior near threshold, the cross-over, and the approach to the saturation value at large p. For definitions 
of symbols and parameter values we refer to the previous figure. 



large p to the mean field result, where density fluctuations are small. 

The corresponding threshold in three dimensional continuum percolation is pp — 0.65296. This 
value is not inconsistent with the low density extrapolation of the three-dimensional heat diffusivity 
in Fig.l, but simulation data for the three-dimensional random solid, sufficiently close to the 
conductivity threshold, are lacking in the simulations of Ref. |24L |25| . 
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Figure 4: Temperature profiles at times of measurement. For p > 3 a steady linear profile is reached, but 
this is not the case for p < 1.6. 



CC 



0.01 



S--..P ' ' ' 


1 ' 


1 ~ 


'■••.□ 




■ 


'■■□ 




. 


'••p 












-D. 


■o 


■ 


■ 


'■Q 

■■a 


■ 








■ 


'A 


■ 


A /n ... 




: 


A/p 


w. 


: N= 10^* 


+ 


"■^- ^.++*- 


■ N=5x10'* 


A 


■■a ^ A AAA' 


■ N= 10^ 


n 


■■■•□aan . 


N=2x10^ 


• 
1 


■"<■■ 



10 



100 



n r^ 



Figure 5: The plot shows that the dominant correction to 1 — Rd{p) behaves in 2D like A2/P with A2 — 1.2. 
The data refer to the deterministic case (vanishing Langevin force) with the Heaviside weight function. 



To further analyze the observed density dependence we show in Fig. 5 the function 1 — R2{p) on 
a log-log plot. The plot strongly suggests that the subleading correction to the heat conduction at 
large p has the form, R2{p) — 1 — A/p with A ~ 1.2, but an analytic calculation of the quantity is 
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effects. The simulations on the top curve (+), middle curve (A) and bottom curve (D) involve resp 
N = 10*, 5 X 10*, 10^ particles. If N increases from 10* to 10^ at fixed p (say p = 50) the system size 
increases from L/vc = \/N/p ~ 14 to 44, and the finite size effects decrease. If A'' and p increase 
by the same factor, the finite size effects remain of the same order, and if the density increases 
from p = 50 to p = 120 at fixed N (say N = 10*), the system size decreases from L/rc — 14 to 
9, and the corresponding finite size effects increase in Fig. 5. The small three-dimensional systems 
seem to be dominated by finite size effects, which appear already before the subleading term A/ p 
becomes dominant. 

It is also of interest to illustrate another effect on the heat diffusivity of the fluctuations. The 
expected coverage with 'black circles' (see Section I) or the black surface fraction at reduced density 
p is (p{p) = 1 — exp[— |;7rp], and the white surface or free volume fraction is exp[— -jTryo] |27l |. Then at 
p — {1.43629; 20; 100} there are in each black circle on average jirp — 1 = {0.13; 14.7; 77.5} excess 
particles, and the white free surface fraction, 1 — (p{p) = exp[— -jTrp] — {0.32; 2 x 10""^; 8 x 10"'^'^}, 
becomes extremely small with increasing p, whereas the corresponding simulated value of the heat 
diffusivity is still a sizeable fraction, 1 — R2{p) — {100%; 10%; 3%}, below its saturation value. 



IV. GENERALIZED HYDRODYNAMICS 

In this section we further explore the analogy between fluids and statistically disordered solids, 
using kinetic theory. Generalized hydrodynamics describes the decay of small spatial fluctuations 
in the conserved local densities at different wavelength, A^ = in/k. Their decay rates depend 
strongly on how the probing wavelength (size of colloidal particles, polymers, pores....) compares 
to the range of the DPD forces. These decay rates can be expressed in terms of fc— dependent 
transport coefficients. The wavelengths cover both the standard hydrodynamic regime as well as 
the mesoscopic regime. 

A classical method J47| to analyze the full hydrodynamic regime is to determine the eigenval- 
ues (decay rates) of the Fourier modes of the linearized Boltzmann equation, and identify the 
/c— dependent transport coefficients from the decay rates. This method is particularly useful if 
there exist intermediate length scales in the problem |4J|, as is the case here. 

The Fourier modes of spatial fiuctuations decay hke exp[— ((fc)t], and may be divided into soft 
(slowly decaying) hydrodynamic modes, and hard (rapidly decaying) kinetic modes. The former 
class has a vanishing decay rate, C(fc ^ 0) oc A;^, in the long wave length limit, and corresponds to 
a conserved density. The latter class consists of hard kinetic modes with a decay rate C,{k — > 0) = 
constant. Generalized hydrodynamics concerns the study of soft modes of fiuctuations of locally 
conserved densities at all wave numbers fc. 

It is the goal of this section to calculate the dispersion relation of the relaxation rate C,{k) of 
the soft heat mode in our model, which corresponds to the locally conserved energy density. The 
behavior of this mode will depend on how the wavelength of the disturbance compares with the 
relevant length scales in the system. 

The relevant length scales in the heat conduction model are the macroscopic system size L, the 
inter-particle distance n"^/'', and the range of the interaction forces r^. The ratio L/rc controls 
the finite size effects, and the reduced density p = nrf is the only dimensionless parameter that 
controls the dynamics of the problem. 

The basic distance to determine whether a perturbation of wavelength Afc decays according to 
standard hydrodynamics with constant transport coefficient is the range Tc, i.e. for Afc ^ re the 
decay of the heat mode is C,{k) ~ k^Dx, where Dt is the standard heat diffusivity. In general 
however, it decays as, 

C(fc) = DT{k)k^ (IV.l) 

with a fc— dependent heat diffusivity Dxik) that approaches the constant transport coefficient Dt 
as k —^ 0. As soon as the wavelength A^ is comparable to r^ the transport coefficient Dxik) 
becomes fc-dependent. This range of excitations is called generalized hydrodynamics. The method 
to study generalized hydrodynamics in the heat conduction model is essentially the same as the 
one followed in Ref. |4J| for a lattice gas cellular automaton model of the Van der Waals equation, 
or in Ref. [32| for standard DPD without energy conservation. 

To set up the kinetic theory we start with the linearized Boltzmann equation. To derive it we 
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So we start with the first equation of the BBGKY hierarchy of the one-, two-, • • • n— particle 
reduced distribution functions, f{xi,t), /(^^(xi,a;2, t), • • • with phase Xi = {ri,ei}. This is done by 
integrating the A'^— particle Fokker Planck or Liouville equation, dtp = Lp with L in l|II.15|l . over 
the phases X2,X3,- ■ ■ , xn- The result is, 

dtfixut) =Jdx2T{12)f^^\xi,X2,t) 

= akBXoJdr2de2di2w{n2)[Ti - T2 + kBTiT2di2]f^H^u^i,ri - R,e2,t) , (IV.2) 



where T(12) is the two-particle Fokker Planck operator, defined through L = J2i<j Tiij) in l|II.15l 



Note that [• • • ] in L of Eas. ljII.15|) have been replaced by [• • • ] in Ea. ljIV.2|) . This implies that the 
correction of relative 0{l/a) to (Ti — T2){1 — 1/a) has been neglected for consistency. Ea. HIV.2p is 
not a closed equation since the time evolution of the one particle distribution function is expressed 
in terms of the pair distribution function. In DPD models with their softly repulsive interactions it 
is in general a reasonable approximation to assume molecular chaos, which expresses the statistical 
independence of the energy fluctuations in different particles [lOj , i.e. 

/(2)(r,e,r',e',i)^/(r,6,i)/(r',e',t). (IV.3) 

The Fokker Planck Boltzmann (FPB) equation for the single particle distribution f{xi,t) is then 
obtained by combining ljIV.2|l and (|IV.3|1 . 

We are specifically interested in studyin g th e decay rates of the Fourier mode, here the heat 
mode, as a function of the wave number k |47l |. So, we study the decay of small deviations from 
thermal equilibrium. This can be done by linearizing the FPB equation around the equilibrium 
distribution function (|II.25|I . /o(x) = •mljQ{e), i.e. 

/(xi, i) = nM^i) [1 + H(^,,t)] , (IV.4) 

where iJ(xi,i) — Hi is a small quantity. This yields, 

dtM^i)Hi ~ n J d^2T{12)M^i)M^2) (1 +Pi2)ffi, (IV.5) 

where the permutation operator, 7^12, interchanges the labels of the two particles i.e. Vijhi — hj. 
Higher terms than first order in H{xi,t) have been neglected. We are interested in the Fourier 
modes of (|IV.5|) . defined as, 

i7(x, t) = e^'^('=)*+*''-'-/i(k, e), (IV.6) 

where k is the wave vector of the Fourier mode and C(fc) its decay rate. The allowed wave numbers, 
ka = iinia/ L with a = x,y,- ■ ■ , and Ua = 0, ±1, ±2, • • • are determined by the periodic boundary 
conditions. Substitution of (jIV.6p into ljIV.5|) yields the eigenvalue equation, 

[ak)+A{k)]Me)h = 0, (IV.7) 

where the operator A(k) is defined as 

A(k)Vo(ei)/i(k, ei)^n f d^2T {l2)M^i)M^2) (l + er'''-'--Vi2) /i(k, ei). (IV.8) 



As a preparation to solve ljIV.8|) we simplify the above expression, by using the relation, 

T{l2)M^i)M^2)B{xiX2) = akl\ow{ri2)di2TiT2M^i)M^2)di2B{xiX2), (IV.9) 

where B{xi,X2) is an arbitrary function of the phases. This simplification combined with the 
relations: e^ — aksTi and wq = pXo gives for the collision operator (|1V.8|I . 



A(k)V'o(ei)/ii = (wo/«) / d€2di2eie2M^i)M^2)di2 [1 + W{k)Vi2] hi, (IV.IO) 

where the Fourier transform of wlR) is. 
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What remains to be done is to solve the eigenvalue equation (|IV.7|I . It is a simple matter to verify 
that h{k — > 0, e) = 1 + a — /3e is an exact eigenfunction of A(k), and substitution in (JIV.IOII yields 
the eigenvalue, 



i(i?^)t„fc2 + 0(fc4). For large k the term 



ak)=oJa[l-W{k)] 

For small k, the k— expansion of W{k) gives W{k) = 1 

W{k) -^ due to the rapid oscillations of exp[— ik • R]. So its limiting behavior is, 

' Drk'^ - Brk* -\ for fcr^ < 1 

u>o for krc ^ 1 



(IV. 12) 



C(fc) = DTik)k^ = 



(IV. 13) 



where Dt is the standard heat diffusivity and Bt the so-called super-Burnett coefficient. Here Dt 
coincides with the mean field result Doo in HIII.10|) - (|IV.12|1 . and 



Ijt — 



LOq 



8d{d + 2) 



(i?4 



(IV. 14) 



with ujQ = pXq. The form (a -|- 1 — /3e) of the eigen mode confirms that this mode is indeed the 
soft mode of interest, corresponding to the conserved energy. 




Figure 6: (a) Decay rate of the 3-D heat mode ^(fc) — DT{k)k in units of cuo and (b) heat diffusivity 
Drik) = Dt — k^Br + ■ ■ ■ in units of ojo^c, plotted versus krc- Solid lines in (a) and (b) are calculated 
with the Lucy function in liii.llll . and a = 100. 



At large k the heat mode becomes a hard kinetic mode with a constant decay rate, (^{k —^ cxd) 
^constant. This behavior is shown in Fig. 6, where the relaxation rate C(fc) and the generalized 
heat diffusivity Drik) in (|1V.12|I are plotted versus kvc for the three-dimensional case. The Fourier 
transforms W{k) can be calculated by using the following formulas. 



W{k) = 



f^ dxxw{x)Jo{qx) 

J„ dxxw{x) 

r, dxxw{x)sin.{qx) 

q /q dxx'^w{x) 



{d = 2) 
id = 3) 



(IV. 15) 



where q — krc and Jo{x) is the zeroth order Bessel function. Note that W{k) is a non-decreasing 
function of k for the Lucy function, and an oscillating one for the Heaviside function. 

It is also worthwhile noting that the calculations of Dt in this section and those in Section IIIB 
give identical results, although the structure of the calculations is very different. In Section IIIB the 
heat flux is calculated in a state of local equilibrium, which is fully factorized because conservative 
forces are absent. In this section on the other hand we have derived a Boltzmann equation for 
the single particle distribution function H{xi, t), based on the molecular chaos assumption, and we 
solved the eigenvalue equation for the heat mode to obtain the eigenvalue or decay rate ({k) ~ k^DT 
at long wavelength. As the molecular chaos assumption also neglects the spatial correlations 
between colliding particles, the kinetic theory results are also mean free results. As discussed in 
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V. CONCLUSIONS AND PROSPECTIVES 

One of the most interesting prospects of the present paper is that the DPD sohd and the Lorentz 
gas are relatively simple many-particle systems that can be used to develop kinetic equations for 
classical fluids that go beyond the mean field Boltzmann equation. Moreover, they are complemen- 
tary to one another in so far as the mechanisms of transport are concerned. In the Lorentz gas it 
is kinetic transport. In the DPD solid the mechanism is collisional transfer. In classical fluids both 
mechanisms are present, and the former is dominant at low and moderate densities; the latter is 
dominant in dense liquids. 

In the Lorentz gas many-particle resummation techniques have been developed that account for 
the collective effects of ring collisions, i.e. sequences of dynamically correlated binary collisions, 
that led to the log-p dependence of transport coefficient 1] and to the power law long time tails 
in the velocity autocorrelation function _2J . 

In the present DPD model we have firmly established by elaborate computer simulations, that 
mean field theory gives essentially exact results for the transport coefficients at very high densities, 
and that the faster-than-linear decrease of the heat diffusivity, Dt{p) = I?oo(p){l — A/ p -!-•••}, 
is caused by local density fluctuations. This was done by comparing on- and off-lattice spatial 
configurations of particles (see Fig.l). 

Consequently, it is to be expected that applications of Effective Medium Theory [43, |43l , which 
is equivalent to the self-consistent ring- kinetic equation 50], and its extensions to classical fluids, 
would provide systematic methods for calculating transport properties, starting at the high density 
side of the density spectrum. 

Furthermore, as the density p decreases the heat diffusivity rapidly decreases to zero at a thresh- 
old density pc, below which the heat conductivity vanishes. The existence of a heat conduction 
threshold pc is explained as a dynamic percolation phenomenon, and identified with bond perco- 
lation on a random proximity network. The dynamics on this discrete network is different from 
diffusion on the continuum percolation structures, although the geometrical connectivity properties 
of the discrete and continuous percolating cluster are the same. The threshold is identified with 
the two-dimensional percolation threshold pc — 1.43629 |23| of continuum percolation of overlap- 
ping spheres, and its value agrees well with the known conductivity threshold pc in Fig. 3a. So the 
simulation values for the heat conductivity agree both at high and low density with our theoret- 
ical analysis of the heat conductivity. Our older three-dimensional simulation data for the heat 
conductivity of Ref. |2J|, shown in Fig.l, are not inconsistent with the existence of a conductivity 
threshold pc at the 3-D percolation threshold 0.65296 z7], but simulation data are lacking close to 
the percolation point, and presumably show strong finite size effects, caused by the small systems 
used in the simulations. 

We have further extended the kinetic theory to the regime of generalized hydrodynamics by 
studying the wave number dependent decay rates of the Fourier modes of the temperature fluc- 
tuations. The decay rate C(fc) = k'^Dxik) of these modes depends strongly on how the probing 
wavelength (size of colloidal particles, polymers, pores....) compares to the range of the DPD 
forces. So, there are several possibilities for applying the generalized hydrodynamics results, apart 
from the applications, already discussed in Section IV, in mode coupling theories, and in analyzing 
finite size effects where the discreteness of the allowed k— values has to be taken into account. In 
the limit of long wave lengths the heat diffusivity Drik -^ 0) reaches the constant value given by 
the standard Chapman Enskog theory. When the wave length of the perturbation is of the same 
order of magnitude as the range re of the forces, the heat diffusivity, predicted by the generalized 
hydrodynamics, decreases significantly below its long wave length value. Here we also mention 
the application of our fc— dependent transport coefficients in smoothed DPD [51|. The goal of such 
methods is to discretize macroscopic nonlinear partial differential equations - here Fourier's law for 
heat diffusivity - and to solve them with molecular dynamics codes (see Ref. [5^ ) . The finite size 
effects, discussed in Section IIIB, are in smoothed DPD, as well as in the related Smooth Particle 
Hydrodynamics ^JIj measures for controlling the discretization errors. The authors of Ref. |5]| 
have measured and analyzed the decay rates (^(k) of sinusoidal temperature profiles in our heat 
conduction model, from which the values of I?T(k) are extracted, and compared with our results 
for DT{k), as presented in Section IV. 
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Appendix A: H-THEOREM AND EQUILIBRIUM STATE 

In this appendix we prove an Ti theorem, and analyze the equilibrium distribution. We show 
that the function 7i is a Lyapunov functional with dtH < 0, and investigate the implications of 
this result for the equilibrium solution of the Fokker Planck equation. 

We consider the following functional of the A^-particle distribution function p(X), 



n[p] = f dx 



Inp(X)- 



5(X) 



P(X) 



(A.l) 



^N 



where 5'(X) = J2i s{ei) and s(ei) is the one particle entropy function, and —Ti. is the total entropy 
of the A^— particle system. Similar results have been obtained in Refs. pH 124. l53| . 

The time derivative of the functional in (jA.ip is given through the Fokker Planck equation, 



dtTi^JdX Inp(X)-^ 9tp(X) == / dXp(X)Lt Inp(X) 
= -/dXp(X)E,<,^.,B.,, 
where we have performed partial integrations, and introduced the symbols. 



5(X) 



(A.2) 



A, 



9. 



ln/?(X) 



5(X) 



d,j In p 



Bi 



The strategy is to make the factor {• ■ • } in Bij equal to Aij, yielding, 

dtH = -JdXp{X) Y, ia2fe-)(Ay)2 < 0, 



(A.3) 



(A.4) 



Kj 



which guarantees that dtTi in ljA.2|) is non- decreasing. This is the desired 7Y-theorem. There are 
two possibilities to realize this. The first one is to choose, 

a^{ij)^2kBX{t]mT, and d,ja^iij) = (V{i,j}). (A.5) 

These conditions are referred to as Detailed Balance conditions. A solution of the last equation is. 



a^{ij) = 2kB\{i])T,T, = K{ij) = 2kBw{r,j)FiT, + T,). 



(A.6) 



Here w{r) is the range function defined in (|II.3p . which implies Wsir) — y/w{r), and F(x) is some 
positive function of x, e.g. F{x) — kqx^"^ in = 0, 1, 2, ...). This is the solution, used in |24L 125115^ . 
With the help of (|II.3(I and l)A.6(l the temperature relaxation equation would take the form, 



dTi 
IT 



Kq 



akBTiT2 



(Ti+T2)^"(T2-ri 



(A.7) 



with n = 0, 1, 2... Although mathematically acceptable this temperature relaxation equation is not 
in agreement with irreversible thermodynamics. 
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Next we discuss the second set of solutions. To do so we drop the second requirement in ljA.5|) 
and write the term {...} in (jA.3|l as, 

{••■} = ^ (^ - ^) +%lnmr,A(jj>] 
1 / 1\ / 1 1 



As explained below (|II.6|I a is a large number, measuring the number of internal degrees of freedom 
of a DPD particle, and (1 — 1/a) should be replaced by 1 for consistency to leading order in \/a. 
Now the expression {...} in l|A.8|) can be made equal to Aij in (|A.2p by choosing \{ij) independent 
of the internal energies of the interacting particles, and in agreement with the laws of irreversible 
thermodynamics 1)11. 12|l and l)II.13() . Then the second set of Detailed Balance conditions becomes, 

a^{ij) — 2kB^{ij)TiTj and A(ij) = akBXow{rij). (A. 9) 

We also quote for later reference that both sets of Detailed Balance conditions lA.Sp and (|A.9p 
guarantee the commutation relation, 

a'(*j)% = d,ja\ij) + 0{l/a). (A.IO) 

The 7i-function in l|A.4|) keeps decreasing and reaches a minimum, if and only if the partial 
differential equations, Aij — 0, are satisfied for all {i,j}- The solution of these differential equa- 
tions determines the equilibrium distribution. As the particles are point particles, the iV— particle 
distribution factorizes. 



Pe,(X) = inVo(Q). (A.ll) 



z 

Combination of these differential equations, with (|A.lip yields for all possible pairs (ij), 

dlnipo 1 QlnV'o 1 



de., ksTi dtj ksTj 



-/3, (A.12) 



where /3 is a constant with dimensions of an inverse energy. By using the relation, ds{ei)/dei — l/T^, 
integration of l|A.12|l yields 

V'o(ej) = -^7g\9i^i) exp[-/3ej] = —— exp[kg^s{e^) - f3£i\, (A.13) 

where the factor g(ei) = exp[fc^ •s(ci)] '^ ^f i^ the degeneracy factor oi the number of internal states 
having energy e^. In a mesoscopic picture g{ei) is the number of internal states of the mesoscopic 
particle i having energy e^. The normalization factor z(/?) is, 

/•OO 

z(/3)=/ decxp[kg^s{e)- f3e]. (A.14) 

Jo 

This factor corresponds to the partition function of a single mesoscopic particle. 

For the DPD solid, where the entropy is given in l)II.5|) . the one particle equilibrium distribution 
function becomes, 

^o(^) = T^^^ (/?^)"exp[-/36], (A.15) 

i (a + Ij 

where r(a;) is the gamma function. For later use we also quote the moments of ipo, 

(,") ^JdeMe)e- = ^^^^^^^ (A.16) 

in particular 
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The parameter (3 is related to the total energy of the system E through the relation, 

E 91nz(/3) ja + l) 

The last equality has been calculated by considering the entropy function ljll.5|) . Furthermore, the 
following average can be calculated 



1 



W) r *^ ''''^ ^^^'"^'^ " ^'^ " ^""^^ ^^'^^^ 



Note that this relation is valid for a general entropy function. It allows one to define the macroscopic 
temperature T as T^^ = (T,^^^) where /? = l/fc^T is the inverse macroscopic temperature in 
thermal equilibrium. For the special choice of IjII.Sp . it is also interesting to point out that, 
(Tj) = [{a + l)/a] T ^ T, this is the relation between the macroscopic temperature and the 
average temperature. As discussed below (jll.611 . a scales like the number of internal degrees of 
freedom and therefore a 3> 1. 

The discussion above deals with fluctuations in the single particle energies, calculated in the 
canonical ensemble. In [2J, |2a| it was shown that such fluctuations, calculated in the micro- 
canonical ensemble, give the same results, provided a is large. This condition is always satisfied 
as DPD particles are mesoscopic objects. 
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